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Abstract 

In this paper, we combine concepts of the generalized multiscale finite element method and mode 
decomposition methods to construct a robust local-global approach for model reduction of flows 
in high-contrast porous media. This is achieved by implementing proper orthogonal decomposi- 
tion (POD) and dynamic mode decomposition (DMD) techniques on a coarse grid. The resulting 
reduced-order approach enables a significant reduction in the flow problem size while accurately 
capturing the behavior of fully resolved solutions. We consider a variety of high-contrast coeffi- 
cients and present the corresponding numerical results to illustrate the effectiveness of the proposed 
technique. This paper is a continuation of the first part [ll where we examine the applicability of 
POD and DMD to derive simplified and reliable representations of flows in high-contrast porous 
media. In the current paper, we discuss how these global model reduction approaches can be com- 
bined with local techniques to speed-up the simulations. The speed-up is due to inexpensive, while 
sufficiently accurate, computations of global snapshots. 

Keywords: Model reduction, generalized multiscale finite element method, heterogeneous porous 
media, dynamic mode decomposition, proper orthogonal decomposition. 

1. Introduction 

Many relevant engineering and scientific applications in porous media problems consist of cou- 
pled processes in highly heterogeneous media. For instance, media permeability can vary over sev- 
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eral orders of magnitude due to high-conductivity fractures and/or low-conductivity shale layers. 
Because of these variations, the resulting large number of degrees of freedom and associated com- 
putational costs could limit the capability to perform a sensitivity analysis or conduct uncertainty 
quantification studies which require solving the forward problem many times. This presents the 
need to develop simplified models that significantly reduce the number of degrees of freedom by 
neglecting irrelevant details of the involved physics in order to remain computationally tractable. 

Local multiscale methods. Multiscale solution techniques represent a class of methods that 
capture the effects of small scales on a coarse grid [zl-ZI- In this paper we follow the framework 
of the multiscale finite element method (MsFEM), where precomputed multiscale basis functions 
are used to span a coarse-grid solution space. As fine scale is embedded into the basis functions, 
we can recover relevant fine scale information from the multiscale solution representation. In re- 
cent years, MsFEM has been extended to allow for the systematic enrichment of coarse solution 
spaces in order to converge to the fine-grid solution (see e.g., llsl-lllll). In addition, the enriched 
coarse spaces have been shown to be effective preconditioners in two-level domain decomposition 
iterative procedures [llOl lllll . The flexibility associated with the construction of enriched coarse 
spaces makes the method an ideal solution technique for a wide variety of applications. In partic- 
ular, the coarse space sizes and associated errors may be carefully calibrated in order to achieve 
appropriate levels of solution accuracy and computational efficiency. In situations where a higher 
level of accuracy is desired, additional basis functions may be incorporated in the construction 
of the coarse space. However, when computational efficiency is a main consideration, less basis 
functions may be used for a significant reduction in the coarse space dimension. In this paper we 
apply the enriched coarse space construction from the generalized multiscale finite element method 
(GMsFEM) as an effective tool for local model reduction [|8l- lll|] . 

Global multiscale methods. Two common techniques have been widely used for global model 
reduction, namely dynamic mode decomposition (DMD) and proper orthogonal decomposition 
(POD). Both of them are based on processing information from a sequence of snapshots (or in- 
stantaneous solutions) to identify a low-dimensional set of basis functions. These functions are 
then used to derive a low-dimensional dynamical system that is typically obtained by Galerkin 
projection [Il2l-ll6ll. Proper orthogonal decomposition constitutes a powerful mode decomposition 
techniqtie for extracting the most energetic structures from a linear or nonlinear dynamical pro- 
cess 1112 - 
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23h . Dynamic mode decomposition has been recently proposed by Schmid [|24ll . 



In comparison with POD, this technique is intended to accurately extract the coherent and dynam- 
ically relevant structures rather than selecting the dominant modes that capture most of the flow's 
energy. DMD enables the computation, from simulation and empirical data, of the eigenvalues and 
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eigenvectors of a linear model that best represents the underlying dynamics, even if those dynamics 
are produced by a nonlinear process. One important feature of this method is its ability to extract 
dynamic information from flow fields without depending on the availability of a model, but rather 
is based on a sequence of snapshots. As such, this technique has been successfully applied to the 



analysis of experimental [|25N30n and numerical [|24l 13 11-13411 flow field data and has shown a great 
capability to capture the relevant associated dynamics and help in the characterization of relevant 
physical mechanisms. 

This paper. In Part I we followed a global approach to derive reduced-order models for 
flows in highly-heterogeneous porous media. This is performed by applying POD and DMD on 
standard finite element solutions. The DMD-based approach showed a better predictive capability 
due to its ability to accurately extract the information relevant for long-term dynamics, in partic- 
ular, the slowly-decaying eigenmodes corresponding to largest eigenvalues. Our contribution in 
this work is to present a local-global model reduction framework for multiscale problems. Our 
approach is based on applying the aforementioned mode decomposition methods to coarse-scale 
problems in order to achieve a significant reduction in the system size while preserving the main 
flow features. In particular, we combine the concepts of GMsFEM with DMD and/or POD to de- 
velop a robust local-global approach for solving high-contrast, time-dependent parabolic problems. 
The resulting reduced-order method is shown to accurately capture the behavior of fully resolved 
solutions for a variety of high-contrast coefficients. Because the accuracy of both local and global 
approaches depends on the number of local and global modes, when selecting local modes, our 
motivation is to obtain inexpensive global snapshots while preserving the accuracy of the global 
approach. In the paper, we achieve this based on a limited number of runs. More extensive apos- 
teriori error estimates can guide a correct choice for local and global modes. We use a variety of 
numerical results to illustrate the effectiveness of the proposed technique. 

The remainder of the paper is organized as follows. In Sect. |2]we describe the problem set- 
ting, after which we describe the local multiscale approach for model reduction in Sect. |3] Global 
model reduction techniques are then introduced and described in Sect. IH In the same section, 
we also present the main steps of the local-global approach that is based on implementing mode 
decomposition methods on the coarse-scale problem. In Sect. |5] we present a variety of numeri- 
cal results to complement the proposed method along with some concluding remarks in the final 
section. 
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2. Model Problem 

In this paper we consider a time-dependent, single-phase porous media flow governed by the 
following parabolic partial differential equation 

— - V ■ («:(x)Vw) = /(x) in n, (1) 

where u denotes the pressure, is a bounded domain, / is a forcing term, and k{x) is a positive- 
definite scalar function. The coefficient k represents the ratio of the permeability over the fluid 
viscosity and is considered to be a highly-heterogeneous field with high contrast (i.e., there are 
large variations in the permeability). The structure of k, is an important factor within this paper, 
and numerous examples will be offered in subsequent sections. We note that Q will be solved 
along with specified boundary and initial conditions. 



3. Local Multiscale Model Reduction 

In this section we describe a systematic coarse grid solution technique that may be used as a 
reduced-order alternative to a standard fine grid approach (such as finite element discretization). 
The solution procedure is built within the framework of the Generalized Multiscale Finite Element 
Method (GMsFEM), where the solution is sought in a space of precomputed multiscale basis func- 
tions. A notable distinction of GMsFEM is that the associated coarse space may be systematically 



enriched to achieve a desired level of numerical accuracy [iSNlOll. In particular, the dimension of 
the coarse solution space may be reliably chosen based on the nature of the problem (e.g., the 
structure of k{x)). The predictable accuracy of the method, combined with the inherent gain in 
efficiency, make GMsFEM a tractable approach for solving the model problem robustly. 

3.1. Fine grid approach 

In order to outline the fine grid solution technique, we first introduce a coarse discretization 
and assume that T'^ is a refinement of T^. To solve © using the finite element method (FEM) 
we search for Uh{t) G = span{(/)j}^^^, where 0^ are the standard bilinear finite element basis 
functions defined on T'\ and Nf denotes the number of nodes on the fine grid. After multiplying 
the equations by test functions and integrating over the domain f2, we obtain the following set of 
ordinary differential equations corresponding to the model equation given in ([I]): 

+ AU = F, (2) 

dt 
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where U = [ui(t)] denotes the time-dependent nodal solution values, M is the mass matrix given 



by M = [rriij] = / A is the stiffness matrix given by A = [aij] = / kV^i ■ V0j, and F is 

Jn Jn 

the forcing vector given by F = [/j] = / Using the backward Euler, implicit scheme for the 

Jn 

time marching process yields 

U'^+i = (M + At A) "^M + (M + AtA)"^ At F, (3) 

where n denotes the time stepping index, and At is the time step. The fine discretization yields 
large matrices of size NfxNj which may become prohibitively expensive to numerically handle. 

Remark 1. Part / ji/ combines the standard FEM solution technique in ([3]) with methods for 
global model reduction. However, the focal point of the present work is to systematically reduce 
the dimension of the system that results from the fully-resolved (fine grid) solution. As this stage of 
dimension reduction hinges on the localized construction of basis functions that span the coarse 
solution space, we refer to it as a local approach for model reduction. By combining methods 
for local and global model reduction we offer a robust framework that efficiently and accurately 
captures the dominant features of the dynamical system for long-term forecasting. 

3.2. Local multiscale approach 

We use {yi}^J^^ (where << Nj) to denote the vertices of the coarse mesh and define 
the neighborhood coi of the node yi by 



|J{i^, er^; y^eK,}. (4) 



See Fig. [T]for an illustration of a coarse neighborhood cui. Using the coarse mesh we start with 
an initial coarse space VJj"'''^' = spanjxj}^]^, where the Xi standard multiscale finite element 
partition of unity functions satisfying 

- V ■ (k(x)Vx*) =0 Ke Ui (5) 
Xi = 9i on dK, 

For simplicity, for all K E cOi, and gi is assumed to be linear. 

A solution computed within V^™'"*' is the standard multiscale finite element (MsFEM) solu- 
tion. However, we may systematically enrich the initial coarse space using generalized approaches 



for multiscale model reduction [|9-[ll|. In particular, we may multiply solutions resulting from 
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Figure 1 : Illustration of a coarse neighborhood 

localized eigenvalue problems to the initial partition of unity to enrich the coarse space. The con- 
struction of the enriched space yields that the solution error decreases with respect to the localized 
eigenvalue behavior. We refer the interested reader to [9l for rigorous error estimates. 

To enrich the initial coarse space, we construct the pointwise energy of the original basis func- 
tions by setting 

|2 



(6) 



i=l 



where H denotes the coarse mesh size. Once k is available, we solve homogeneous Neumann 
eigenvalue problems of the form 



(7) 



on each coarse block neighborhood coi. We denote the eigenvalues and eigenvectors of (|7]) by 
{A^'} and {'ip'^'}, respectively. Since we consider a zero Neumann problem, we note that the first 
eigenpair is X^' = and = 1. We order the resulting eigenvalues as 



Af < A2' < . . . < XT < 



(8) 



The size of the eigenvalues is closely related to the structure of k and, in particular, m inclusions 
and channels in k yields m asymptotically vanishing eigenvalues. It is precisely the eigenvectors 
corresponding to small, asymptotically vanishing eigenvalues that we wish to use for the construc- 
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tion of the coarse space Vq. As such, we define the basis functions 

$M = Xii^r for 1 < 2 < N, and 1 < / < L^, (9) 

where Lj denotes the number of eigenvectors that will be chosen for each node i. With the updated 
basis functions in place, we define the local spectral multiscale space as 



Vo = span{<l>,,i : 1 < i < A^^nd 1 < / < L,}. (10) 

Using a slightly different index notation, we may write Vq = span{$j}^^, where denotes the 
total number of basis functions used in the coarse space construction. 

Through constructing an operator matrix Rq = . . . , $ArJ (where are taken to be the 
nodal values of each basis function defined on the fine grid), we can express the coarse scale 
analogue of © as 

U^i= (Mo + AtAo)"'MoU^ + (Mo + AtAo)-'AtFo, (11) 

where Mq = RqMR^, Aq = RqAR^, and Fq = RqF. To elaborate on the form in (fTTTl . 
we reiterate that we now seek solutions within the spectral multiscale space Vq. A more detailed 
consideration of the resulting coarse scale matrices (using the mass matrix as a specific example) 
yields an expression of the form 

Mo = = / = RoMR^, 

Jn 

since for all p,q eV^ we have p^yiq = / pq. Rq analogously allows for the extension of coarse 

Jn 

scale solutions onto the fine grid. The resulting coarse matrices in (111]) are of size NcXNc, where 
Nc is significantly smaller than Nj. Thus, the coarse system in (fTTI) offers a suitable local model 
reduction of the fine system in (|3]). 



4. Global Model Reduction 



Two mode decomposition methods are used in this study for global model reduction, namely 
Proper Orthogonal Decomposition (POD) and Dynamic Mode Decomposition (DMD). The basic 
principles of POD and DMD are presented in this section. For more detailed description of POD, 



the reader is refereed to [jl2l-ll5l 117^2311 and for DMD DM IZSl InL l35|] . 



The first step for either the POD- or DMD- analysis is to collect a sequence of instantaneous 
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data fields (or snapshots) \j given by 



Vf = {Vi,V2,V3,--- ,v^} (12) 

where the time spacing between two consecutive snapshots in the above sequence is assumed to 
be constant. The POD proceeds by performing a singular value decomposition of the sequence of 
snapshots V^. To this end, we first form the correlation matrix C from the snapshot sequence as: 

C = (Vf )^ Vf , (13) 

and then compute the POD modes 0^*^^ by performing an eigen-analysis of the correlation matrix 
C; that is, 

CWi = al Wi and 0f^^ = — Vf Wi. (14) 

The selection of POD modes is based on an energy ranking of the coherent structures. However, 
the energy may not in all circumstances be the appropriate measure to rank the importance of the 
flow structures and detect the most dynamically -relevant modes III Bsll . 

The DMD method is based on postprocessing the sequence of snapshots \^ to extract the 
dynamic information. It uses the Amoldi approach to relate two consecutive data fields through a 
linear mapping; that is, 

Vi+i = Avi. (15) 

The eigenvalues and eigenvectors of the matrix A characterize the dynamic behavior of the flow. 
In practical situations, the matrix A might be very large or even its exact form may not be given. 
As such, the DMD method enables the approximation of the relevant eigenvalues and eigenvectors 
of the matrix A. To proceed with DMD, we first assume that the last snapshot can be represented 
by a linear combination of the previous snapshots; that is, 

Af-l 

Vat = ^ fliVi + r (16) 

i=l 

or 

VAr = Vf-ia + r (17) 
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where a = {ai, 02, ■ ■ ■ > O'N-i}'^ and r is the residual vector. This yields 



A Vf -1 = = Vf S + r e^_i 



(18) 



where e^^^^ ~ ( ' ' ' ^ 1 j is the (iV — 1) unit vector and the matrix S is of companion type 



defined as: 



/o 

1 



v 



0-2 



1 aN^2 
1 aAT-i J 



(19) 



The unknown matrix S is determined by minimizing the residual r. The minimization problem 
expressed as 



min I V 
s 



(20) 



can be solved using a QR-decomposition of the matrix Vf^~^. Once S is determined, we com- 
pute its eigenvalues and eigenvectors can be computed which result in the DMD modes 0^*^^ as 
follows: 



S = R"^ 



(21) 



^DMD ^ yN-1 (22) 

where 

SXj = XjXj. (23) 



A more robust approach to compute the DMD modes [|24|1 is based on a preprocessing step using 
a singular value decomposition of the data sequence Y^~^ = U S W^. Substituting the SVD 
representation U S into Eq. (fTSi) and multiplying the result by from the left and by W 
from the right, we obtain the following matrix 

A U = W = S. (24) 



9 



The dynamic modes are then computed from the matrix S as follows: 



^r'' = uy„ (25) 



where is the j*-^ eigenvector of S, i.e., S = P-jYj, and U is the matrix collecting the right 
singular vectors of the snapshot sequence Vf^~^. The latter approach is used in the present study. 

5. Numerical Results 

In this section, we analyze the effects of permeability fields that contain inclusions of high and 
low conductivity as shown in Fig. [21 These configurations lead to different types of dynamical 



behavior. Fig. |2(a)| is a field that models high conductivity channels within an otherwise homo- 
geneous domain. The minimum conductivity value for this case is taken to be = 1, and 
the high conductivity channels vary randomly (from channel to channel) with a maximum value 
of Kmax = 1.9x10^. Fig. |2(b)| is a field of similar structure, yet the permeability values are in- 
verted such that the layers (similarly placed) now represent low conductivity regions within a high 
conductivity value. Both permeability fields are described by a 100x100 fine mesh. Using such 
configurations, we solve the problem governed by Eq. (fTTT i on a coarse grid over a time interval 
of T = [0 1] . This time interval is observed to be long enough so that the steady-state solution 
for all relevant cases is achieved. We use a time step of At = 5x 10"^, and solve the model on a 
two-dimensional unit domain f2 = [0, l]x[0, 1]. We assume zero Dirichlet boundary conditions, 
and the initial solution configuration is shown in Fig. |3l 

5.1. Local multiscale model reduction results 

To further motivate the combined local-global solution procedure, we first present results from 
the local multiscale approach within this subsection. As the local-global approach described in 
Sect.|4]does not require solutions that are obtained through a coarse grid approximation (see 111), 
the addition of an effective multiscale model will offer a more efficient solution technique which 
is robust. For these initial comparisons we use a forcing term / = 1, the permeability field from 
Fig. |2(a)[ and recall that the fine solutions are obtained from Eq. ([3]) and the multiscale (coarse) 



solutions are obtained from Eq. (fTTT) . For these examples, the fine mesh yields a system of size 
Nf = 10201, and the local multiscale approach yields a system of size Nc = 850. Thus, we 
are solving a reduced-order system which is an order of magnitude smaller than its fully-resolved 
counterpart. See Fig. |4]for an illustration of fine and reduced-order solutions advancing in time. 
The solution profiles are nearly indistinguishable from one another, formally validating the success 
of the local multiscale approach. For a more rigorous comparison, we also offer the relative L2 
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(a) Case I: high-conductivity channels 



(b) Case II: low-conductivity layers 



Figure 2: Different configurations of the permeability field; high conductivity channels (left), low conductivity layers 
(right). 

error quantities \\u^f{t) — u^''{t)\\ / \\u^^{t)\\ x 100% for a representative amount of time in 



The error values in Fig. |5]typically range from 1.0 — 1.5% for any time value. This range of 
values results from a careful choice of the coarse space dimension that we use in the local model 
reduction. In particular, we choose the dimension of the coarse space such that the errors associated 
with GMsFEM are comparable to those that will result from the respective global model reduction 
technique. In doing so, we avoid the larger computations that are associated with achieving a 
more strict level of accuracy (due to the addition of more basis functions in the coarse space 
construction). The alternative would be to devote unnecessary resources to obtain solutions whose 
errors would be essentially negligible compared to those those resulting from the global model 
reduction. So while it may be possible to compute more accurate solutions within a larger coarse 
space, such an effort would be in excess of the optimal execution of the proposed method. 

5.2. Local-global model reduction results 

The variations of the forcing term / in our numerical examples over Vt is shown in Fig. [6l The 
objective is to investigate the capability of a combination of global model reduction techniques, 
POD and DMD, with local multiscale approach to capture the main flow characteristics and repro- 
duce the flow dynamics response within a certain accuracy and at a reduced computational cost. 
The local-global approach involves the following steps: 

• Using the same solution parameters from above, we record Nt instantaneous solutions (usu- 



Fig.El 
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Figure 3: Initial configuration of the solution field, 
ally referred as snapshots) and collect them in a snapshot matrix as: 

Vf = {Vi,V2,V3,---,VjvJ (26) 

where Nt is the number of snapshots and is the size of the column vectors Vj. A particular 
attention should be drawn when selecting the set of snapshots that will be used to extract the 



modes fll)]. For both cases, we start collecting snapshots from the 10*-^ time step. The 
use of 25 snapshots is observed to be appropriate to reproduce results of Case I (inclusions 
with high conductivity) with acceptable accuracy. However, 100 snapshots were required 
to properly analyze Case II. This is expected since the permeability field that corresponds 
to Case II contains inclusions of low conductivity and subsequently involves longer time 
dynamics. 

We postprocess the snapshot matrix, as described in the previous section, to compute the 
POD and DMD modes and use these modes to approximate the solution field in the coarse 
grid. As such, we assume an expansion in terms of the modes 0,^; that is, we let 

Nr 

Uo{x, t) ^ Uo{x, t) = ^ 0!i{t)(l)'l{x) ill) 
i=l 
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Fine Solution (t = 0.01) Fine Sohition (t = 0.02) Fine Sohition (t = 0.03) 




(a) Fine (Nf = 10201) 



Coarse Solution (t = 0.01) Coarse Solution [t = 0.02) Coarse Solution (f = 0.03) 




XXX 



(b) Coarse (N^ = 850) 
Figure 4: Fine (top) and coarse (bottom) solutions advancing in time 

or in a matrix form 

UJJ ^ = (28) 

where $ = ^ 0^ ■ ■ ■ j and k can refer to POD or DMD. For the sake of comparison, 
we extract and keep the same number of modes for POD and DMD in our simulations. 
We consider the use of the first six modes for all simulations. In fact, the singular value 
decomposition of the snapshot matrix showed that the first six POD modes contain more 
than 99.9% of the total flow energy. 

• To assess the capability of the POD and DMD modes in capturing the dynamics involved 
in the process and enabling good projection subspaces, we evaluate the L2 projection error; 
that is, we project each snapshot onto the POD modes, compute the following inner product 

ai{t,) = (0(=(x),Mo(x,t,)) , or = ($*<l>)-i$*U^o (29) 

where (F, G) = Jq{F G) dVt, and define the relative error as the L2-norm of the difference 
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Multiscale Errors 




Figure 5: Relative errors between fine and multiscale solutions advancing in time 



between the reference and approximate solutions over the reference one; i.e. 

II uo{x,t) - uo{x,t) II2 



m 



Uo{x,t) 



(30) 



A low projection error indicates the ability of modal decomposition techniques to compute 
good projection subspaces, but it does not necessarily yield stable and accurate reduced-order 
models when Galerkin projection is used. 



We use the solution expansion given by Eq. (1271) and project the governing equation of the 
coarse scale problem onto the space formed by the modes to obtain a set of Nr ordinary 
differential equations that constitute a reduced-order model; that is, 



a 



($*Mo$)"'$*Ao$a+ ($*Mo$)"'$*Fo. 



(31) 



Thus, the original problem with Nf degrees of freedom is reduced to a dynamical system 
with Nr dimensions where Nr <^ Nf. 

• We use the operator matrix Rq to downscale the approximate solution and evaluate the flow 
field in the fine scale domain. 

We follow the above steps to derive reduced-order models while considering the different initial 
configuration shown in Fig. Ul This is performed intentionally to test the robustness of the reduced- 
order model with respect to variations in the initial conditions. In Fig. [8l we plot the variations of 
the L2 projection and Galerkin projection errors with time for Cases I and 11. Results are obtained 
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Figure 6: Spatial variations of the forcing / over the domain O. 

from the DMD- and POD-based approaches. We observe small L2 projection errors. However, 
larger errors are reached when projecting the governing equations onto the subspace spanned by 
the modes to obtain a reduced-order model. In particular, we observe jumps to high values during 
the first time steps. This is due to the use of different initial configuration. For both cases of high 
and low conductivity, the Galerkin projection errors decrease and reach relatively small values as 
time evolves and the steady state develops (Case I: 12% for POD and 0.6 % for DMD and Case 
II: 14% for POD and 3 % for DMD). This indicates the suitability of the combined local-global 
approach for model reduction of flows in highly heterogeneous porous media. 

Figs. ID and [To] depict the output fields at t = 1 obtained from the reference fine scale solution 
and those obtained from the local-global approach for both high and low conductivity cases. We 
observe good agreement between the reference solution and that obtained from the hybrid DMD- 
coarse multiscale approach. However, a discrepancy can be seen when comparing the reference 
solution with that obtained from the hybrid POD-coarse multiscale approach. These observations 
show the capability of DMD modes, computed from the first few snapshots, to capture the relevant 
flow dynamics and forecast the flow field with good accuracy for long time periods. 

To investigate further the suitability of POD and DMD modes to model flows in varying and 
highly heterogeneous porous media, we consider Case I; that is, the permeability field shown in 
Fig. |2(a)| and multiply its coefficient by a smooth positive spatial function to obtain 

Ks{x; y; e; /) = y) x {1 + e + sin(27r/x) sm{27i fy)), (32) 
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Figure 7: Initial configuration of the solution field. 



where e = 2 and / = 25. The resuking permeability field is depicted in Fig. \TT\ This analysis is 
motivated by several applications which require solving the forward problem for varying perme- 
ability fields; for instance, when the permeability field is subject to uncertainty or multi -phase flow 
where the permeability is modulated by coarse-grid mobility. To address this issue, we use POD 



and DMD modes generated for the permeability field shown in Fig. |2(a)| and employ the Galerkin 
projection to obtain a reduced-order model which is used to predict the flow field resulting from 
the modified permeability field described by Eq. (|32] ) shown in Fig. [TT] In Fig. [121 we plot the 
temporal variations of the Galerkin projection and L2 projection errors obtained from the POD- and 
DMD-based representations. Low steady-state errors are observed. In particular, smaller values 
are obtained when implementing the dynamic mode decomposition procedure on the coarse- scale 
problem. This shows the robustness of the hybrid DMD-coarse multiscale approach with respect 
to moderate perturbation in the permeability field. 



6. Conclusion 

In this work, we propose a local-global approach for model reduction of high-contrast and 
time-dependent parabolic problems that govern flows in highly-heterogeneous porous media. This 
approach combines the concepts of generalized multiscale finite element method (GMsFEM) and 
proper orthogonal decomposition (POD) and/or dynamic mode decomposition (DMD) techniques. 
We consider different high-contrast coefficients and present numerical results to investigate the 
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(a) case I: high-conductivity channels (b) case II: low-conductivity layers 

Figure 8: Variations of the Galerkin projection error (solid line) and L2 projection error (dashed line) with time for 
different permeability configurations. Results are obtained using POD and DMD modes. 




(a) Fine -Reference solution (Nf = (b) Coarse-DMD -Approximate solu- (c) Coarse-POD -Approximate solu- 
10201) tion (Nr = 6) tion (Nr = 6) 



Figure 9: Output fields at t = 1 (Case I: high-conductivity channels). Comparison between reference solution of the 
fine scale problem with those obtained from the hybrid DMD- and POD-coarse multiscale approaches. 



capability of our proposed approach to accurately capture the behavior of resolved solutions. The 
hybrid DMD-GMsFEM technique shows great potential to reproduce the flow field with good 
accuracy while reducing significantly the size of the original problem. This is due to the systematic 
construction of accurate coarse spaces from GMsFEM along with DMD's ability to extract the 
dynamic information and especially the relevant modes that govern the long-time dynamics. This 
achievement opens the door for large-scale optimization applications where the stochasticity of 
the parameters as well as their uncertainty can be efficiently involved in large-scale and accurate 
simulations. 
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Figure 10: Output fields at f = 1 (Case II; low-conductivity layers). Comparison between reference solution of the 
fine scale problem with those obtained from the hybrid DMD- and POD-coarse multiscale approaches. 
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